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Abstract We use our theory of periodized discrete elasticity to characterize defects in graphene as 
the cores of dislocations or groups of dislocations. Earlier numerical implementations of the theory pre- 
. dieted some of the simpler defect groupings observed in subsequent Transmission Electron Microscope 

■ experiments. Here we derive the more complicated defect groupings of three or four defect pairs from 

C/^ I our theory, show that they correspond to the cores of two pairs of dislocation dipoles and ascertain 

^ • their stability. 
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1 Introduction 



o . 

O . Graphene is a two-dimensional (2D) allotrope of carbon formed by a single layer of graphite that was 

first synthesized in 2004 [T]. Since then, graphene has attracted great attention due to the Dirac-like 
^sq ! spectrum of its charge carriers and the resulting extraordinary electronic properties [2j[3]. This fasci- 

^ ■ nating material has remarkable electronic and mechanical properties [2][3j|4]. The electronic, chemical, 

00 I thermal and mechanical properties of graphene are exceptionally sensitive to lattice imperfections [2] 

[3] and these defects and even the ripples that always cover suspended graphene sheets [5]|6] induce 
OO ■ pseudo-magnetic gauge fields [4 . Thus the study of defects in graphene is crucial and it has generated 

\ important experimental work [7.8.9. 10^. Among observed defects, there are pentagon- heptagon (5-7) 

pairs. Stone- Wales (SW) defects (5-7-7-5 defects) [7 , pentagon-octagon-pentagon (5-8-5) divacancies 
[9], asymmetric vacancies (nonagon-pentagon or 9-5 pairs) and more complicated groupings such as 
5-7-7-5 and 7-5-5-7 adjacent pairs or defects comprising three pentagons, three heptagons and one 
hexagon j7]. In other two dimensional (2D) crystals such as Boron Nitride (hBN) symmetric vacancies 
have been observed [11] . 

■ Real time observation of defect dynamics is possible using Transmission Electron Microscopes 

\ (TEM) corrected for aberration that have single atom resolution 0. Defect dynamics in graphene 
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Fig. 1 (Color online) Structure of a graphene lattice. Neighbors of a given atom A in sublattice 1 (dark red). 
Atoms in sublattice 2 are pink. The primitive vectors a, b and c are also indicated. 



occurs on a time scale of seconds [7 , much longer than sub-picosecond time scales typical of sound 
propagation in a primitive cell. On this long time scale and for unstressed graphene, SW defects are 
unstable: their two 5-7 pairs glide towards each other and annihilate, and the same occurs to defects 
comprising three pentagons, three heptagons and one hexagon, whereas 5-7-7-5 and 7-5-5-7 adjacent 
pairs remain stable [7|. In stressed graphene oxide samples, SW defects split into their component 
5-7 pairs which then move apart [10 . While most theoretical studies on the influence of defects in 
electronic properties assume a given defect configuration and then proceed to analyze its effects [3,4 , 
it is important to predict defect stability and evolution. 

In recent work, we have explained the observed long time defect dynamics in graphene by consid- 
ering defects as the core of edge dislocations or dislocation dipoles in a planar 2D hexagonal lattice 
p. 2 pi 3' . Our theory is a top-down approach whose starting point is linear elasticity. We discretize con- 
tinuum linear elasticity on a hexagonal lattice and replace differences of vector displacements along 
primitive directions by periodic functions thereof which are linear for small differences. Our periodized 
discrete elasticity allows dislocation gliding along primitive directions and it reduces to continuum 
linear elasticity very far from dislocation cores [14 . Introducing a large damping in the resulting equa- 
tions of motion and solving them numerically, we are able to predict the stable cores corresponding to 
a given dislocation configuration. Using this theory, we have predicted the stability of 5-7 defects (that 
are the cores of dislocations) [T2j[T3]. Similarly, a study of dislocation dipoles in unstressed samples 
p!2l[T3| predicts that SW are unstable whereas symmetric vacancies, divacancies and 7-5-5-7 defects 
are stable. In stressed samples, our theory predicts that SW split into two 5-7 pairs that move apart 
[13], as confirmed later by experiments [10]. 

In this paper, our theory is used to explain the evolution of defects involving a pair of 5-7-7-5 and 
7-5-5-7 defects and a defect comprising three heptagons, three pentagons and one hexagon as observed 
by Meyer et al [7|. 

The rest of the paper is as follows. Our theory and its equations of motion are explained in Section 
[2j The stable cores corresponding to the far field of a single edge dislocation and a single dislocation 
dipole are used in Section [3] to illustrate the way defects are constructed numerically. Our results 
are also compared to available experiments in graphene and other 2D crystals. Section H) contains 
the new results on more complex defects comprising two dislocation dipoles. We explain experimental 
observations by Meyer et al [T. The last section is devoted to our conclusions. 



2 Periodized discrete elasticity of planar graphene 

In this paper we consider a planar graphene sample and ignore possible vertical deflections. In the 
continuum limit, in-plane deformations are described by the Navier equations of linear elasticity for 
the two-dimensional (2D) displacement vector (u^v) [15^ • Including a phenomenological damping with 
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coefficient 7 (to be fitted to experiments), we have 

where p2 is the 2D mass density and A and ji are the 2D Lame coefficients. 

The governing equations of our theory are obtained in a three step process [13]: (i) discretize (1) on 
the hexagonal graphene lattice, (ii) rewrite the discretized equations in primitive coordinates, and (iii) 
replace finite differences appearing in the equations by periodic functions thereof in such a way that 
the equations remain invariant if we displace the atoms one step along any of the primitive directions. 
The last step allows dislocation gliding. 



2.1 Step (i): Discrete elasticity 

Next we discretize the equations of motion on a hexagonal lattice using the same notation as in Ref. 
[13 . Let us assign the coordinates {x^y) to the atom A in sublattice 1 (see Figure [T]). The origin of 
coordinates, (0, 0), is also an atom of sublattice 1 at the center of the graphene sheet. The three nearest 
neighbors of A belong to sublattice 2 and their cartesian coordinates are ni, 77,2 and 77,3 below. Its six 
next-nearest neighbors belong to sublattice 1 and their cartesian coordinates are n^, i = 4, . . . , 9: 

fa a \ fa a \ f a 

/ a ay/?>\ f a a>/3\ , . 

^4 = \x--,y — , = \x^-,y — , ne = {x-a,y), 



2'" 2 /' V 2'" 2 



a 



a\/?>\ f a a>/3 



nr = (x + a, ns = [x- -,y^ — - ,719= \x^-,y^ — - . (3) 



In Fig. [11 atoms uq and nj are separated from A by the primitive vector ±a and atoms 714 and ng 
are separated from A by the primitive vector ±b. Instead of choosing the primitive vector ±b, we 
could have selected the primitive direction ±c along which atoms ng, A and lie. Let us define the 
following operators acting on functions of the coordinates (x, y) of node A: 

Tu = [u{ni) - u{A)] + [u{n2) - u{A)] + [u{ns) - u{A)] - {d^u + d^u) ^, (4) 

Hu = [nine) - u{A)] + [u{n7) - u{A)] - {dlu) (5) 

/ 1 a/3 3 \ 

Diu = [u{n4) - u{A)] + [u{n9) - u{A)] - f 4 + ^ d^dyU ^ -dyUj (6) 

D2U = [u{n,) - u{A)] + [u{ns) - u{A)] ^ Q d^u - ^ d^^dyU + ^ d'^u^ a^ (7) 

as the lattice constant a tends to zero. Similar operators can be defined if we replace the point A 
in sublattice 1 by a point belonging to the sublattice 2. Now we replace in ^ and (j2j), Hu/a^^ 
(4T — H)u/a^ and [Di — D2)u/ {y/^a^) instead of d^u^ dyU and dxdyU, respectively, with similar 
substitutions for the derivatives of thereby obtaining the following equations at each point of the 
lattice: 

P2a^d^u + 7 Stii = 4/i Tii + (A + /i) Hu + {Di - D2)v, (8) 

v3 

P2a^d^v -^-fdtv = 4(A + 2/i) - (A + p)Hv + {Di - D2)u. (9) 

V3 



These equations have two characteristics time scales, the time tg = \J P20?' I (A + 2/i) it takes a longi- 
tudinal sound wave to traverse a distance a and the characteristic damping time, = ^a^ /{\ + 2/i). 



4 



L.L. Bonilla, A. Carpio 



Using the known values of the Lame coefficients at the graphite basal plane [16^ Q (that agree with 
calculations [T7] and measurements in graphene [H]), ts ~ 10~^^ s. Our simulations show that it takes 
OAtd a SW to disappear after it is created by irradiation which, compared with the measured time of 
4 s gives td ^ 10 s. On a td time scale, we can ignore inertia in (|8])-(j9j). 



2.2 Step (ii): Nondimensional equations in primitive coordinates 



We now transform ([HI)-® to the nondimensional primitive coordinates ^ v' using u = a{u' + '^V^), 
V = \/?>av' use the nondimensional time scale t' = t/td and ignore inertia. The resulting equations 



are 



dt'u' 



A + 2/i A + 2/i 



H 



H 



- Do 



2T]v' 



3 A + 2/i 



A + 2/i 



(10) 
(11) 



2.3 Step (iii): Periodized discrete elasticity 



The models described by the linear equations (p!Q|) - (pT]) do not allow for the changes of neighbors 
involved in defect motion. One way to achieve these changes is to update neighbors as a defect moves. 
Then (p!Q|) and (pT]) would have the same appearance, but the neighbors rii would be given by (|3|) only 
at the start. At each time step, we keep track of the position of the different atoms and update the 
coordinates of the n^. This is commonly done in Molecular Dynamics, as computations are actually 
carried out with only a certain number of neighbors. Convenient as updating is, its computational cost 
is high and analytical studies thereof are not easy. 

In simple geometries, we can avoid updating by introducing a periodic function of differences in the 
primitive directions that automatically describes link breakup and union associated with defect motion. 
Besides greatly reducing computational cost, the resulting periodized discrete elasticity models allow 
analytical studies of defect depinning [l9[|ll], motion and nucleation [20l[2T]. Another advantage of 
periodized discrete elasticity is that boundary conditions can be controlled efficiently to avoid spurious 
numerical reflections at boundaries. 

To restore crystal periodicity, we replace the linear operators T, Di and D2 in (p!Q|) and (pT]) by 
their periodic versions: 

T^u' = g{u\ni) - u\A)) + g{u\n2) - u\A)) + g{u\ns) - u\A)), 
H^u' = g{u\ne) - u\A)) + g{u\nj) - u\A)), 
Dipu' = g{u\n4) - u\A)) + g{u\n9) - u\A)), 
D2pu' = g{u\n^) - u\A)) + g{u\ns) - u\A)), 

where ^ is a periodic function, with period one, and such that g{x) ^ x as x — > 0. We obtain: 



A + 2/i A + 2/i 



Dip — D2p \ , 



In our tests we have taken ^ to be a periodic piecewise linear continuous function: 



. . fx, —a < X < 

g<.{x) = I + j^, a<x<l-a. 



(12) 

(13) 
(14) 

(15) 



^ At 300 K, Cee = Msd = 440 GPa, C12 = Ago = 180 GPa, Cn = A3D + 2/i3D = 1060 GPa. The 2D 
coefficients are A = Xsod and = fisod, where d = 3.35 A is the distance between graphene planes in 
graphite. Similarly, p2 = pd is found from the bulk density of graphite. 
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Fig. 2 (Color online) To generate defects, we use lengths depicted in the figure as referred to a hexagon cell: 
a and / = a/VS are the lattice constant and the hexagon side, respectively. H — I -\- h, where h — 1/2 is the 
vertical distance between nearest neighbor atoms. 

The parameter a controls defect stability and mobility under applied stress. It should be sufficiently 
large for elementary defects (dislocations, vacancies) to be stable at zero applied stress, and sufficiently 
small for dislocations to move under reasonable applied stress [14 . We use a = 0.4 to account for 
experimentally observed stability properties of the defects. For lower values, the stable defect described 
in section m loses the Stone- Wales component. The periodic function g can be replaced by a different 
type of periodic function to achieve a better fit to available experimental or numerical data. 



3 Stable cores of dislocations and dislocation dipoles 

3.1 Boundary and initial conditions for a single dislocation 

We solve (p!3 |) - (p!l|) . or (p!Q]) - (pT]) with the periodic operators Tp, Hp, Dip and I^2p, using as initial and 
boundary conditions the far field of appropriate dislocations which are the stationary solutions of the 
linear elasticity equations [15 . Since the latter are a good approximation four spacings away from the 
core of SW defects in graphene, and our model equations seamlessly reduce to linear elasticity in the 
far field, we use a relatively small lattice with 18 x 18 spacings (36 x 36 carbon atoms) in our numerical 
simulations [13]. Consider first the case of a single edge dislocation with Burgers vector (a, 0) and 
displacement vector u = {u{x,y),v{x,y)) 



where u = A/[2(A + ji)] is dimensionless; cf. Ref. [15 , pag. 114. (p!6|) has a singularity ex {x^ + V^^~^^^ 
at the origin of coordinates and it satisfies /^((ix • V)u = — (a,0), for any closed curve C encircling 
the origin. Using ([16]), we write u = [u^v) in primitive coordinates, U\l,m) = [u{x — xo^y — yo) — 
v{x - xo,y - yo)/VS\/a, V\l,m) = 2v{x - xo,y - yo)/{aV3), where x = {x' ^ y'/'^)ci, y = V^ay'/2, 
x' = I, y' = m (integers) and (xq, ^o) 7^ (0? 0) to avoid that the singularity in (p!6|) be placed at a lattice 
point. To find defects, we solve the periodized discrete elasticity equations (p!3l) - (fT4l) with the initial 
and boundary conditions: 

u'(/, m; 0) = U'(/, m), and u'(/, m; t) = U'(/, m) + F(m, 0) at lattice boundaries. (17) 

Here F is a dimensionless applied shear stress. For |F| < Fc (Peierls stress), the solution of (p!3l) - (p!4|) 
relaxes to a stable dislocation {u' {I , m) , v' {I , m)) with appropriate far field, which is (p!6|) if F = 0. 

Numerical simulations give us the location of carbon atoms at each time t. We represent atoms by 
spheres of arbitrary size. As a guide to the eye and to visualize defects more easily, we have attached 
fictitious bonds to these spheres [12 ^ .13 ^ . Depending on the location of the singularity (xc^o), there 




(16) 
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a) 



b) 



Fig. 3 (Color online) (a) Symmetric vacancy, (b) Asymmetric vacancy (nonagon- pent agon defect). 

are two possible configurations corresponding to the same edge dislocation in the continuum limit. 
If the singularity is placed between two atoms that form any non-vertical side of a given hexagon, 
the core of the deformed lattice (/ + u'{l^m)^m + v'{l^m)) is a 5-7 (pentagon-heptagon) defect. If the 
singularity is placed in any other location different from a lattice point, the core of the singularity 
forms an octagon having one atom with a dangling bond [T2]|T3]. Stable 5-7 defects are commonly 
observed in experiments [71 [TT| [TQ], whereas adsorbed atoms (not considered in our model) may attach 
to a danghng bond thereby destroying the octagon configuration. 



3.2 Dislocation dipoles 

A dislocation dipole is formed by two dislocations with opposite Burgers vectors, ±a. Depending 
on how we place the origin of coordinates, different dipole configurations result. Let U(x, y) be the 
displacement vector (p!6|) of a single dislocation. We find the dipole cores by selecting as zero stress 
initial and boundary conditions U(x — x'^^y — y'^) — U(x — Xq ^y — y^)^ with different {x^^y^). Let a, 
/ = a/>/3, H = 31/2 and h = 1/2 be the lattice constant, the hexagon side, the vertical distance between 
two nearest neighbor atoms belonging to the same sublattice, and the vertical distance between nearest 
neighbor atoms having different ordinate, respectively; cf. Fig. [2j We get: 

— Vacancies: Xq = —0.25a, y^ = —0.8h-\-H/2 and Xq = —0.25a, y^ = —0.8h. This initial configura- 
tion is the asymmetric vacancy (9-5 defect) of Fig. [3l^b), which evolves to the symmetric vacancy 
of Fig. [3l^a) under overdamped dynamics. 

— Stable divacancy: Xq = —0.25a, ^q" = —O.Sh + H and Xq = —0.25a, y^ = —O.Sh. Fig. SI 

— Stable 7-5-5-7 defect: x+ = -0.25a -|- a, y^ = -O.Sh and x^ = -0.25a, % = -0.8/i + iJ. Fig.[5];a). 
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Fig. 5 (Color onhne) (a) Stable 7-5-5-7 defect, (b) Unstable 5-7-7-5 Stone Wales defect. 
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— Unstable Stone- Wales 5-7-7-5 defect: Xq = — 0.25a + a, = —0.8h and Xq = —0.25a, i/q = —0.8h. 
Fig.[5fb). For F = 0, this initial configuration corresponds to two dislocations with opposite Burgers 
vectors that share the same glide line, and it evolves to the undisturbed lattice when the dislocations 
move towards each other and annihilate. 



3.3 Comparison with results of experiments 

Carbon atoms and defects in graphene sheets are visualized by operating at low voltage (< 80 kV, 
to avoid irradiation damage to the sampl^J) a transmission electron aberration-corrected microscope 
(TEAM) with appropriate optics [7 . This microscope is capable of sub- Angstrom resolution even at 
80 kV and can produce real time images of carbon atoms on a scale of seconds: each frame averages 
Is of exposure and the frames themselves are 4 s apart [71122]. The images obtained in experiments 
can be used to determine the time evolution of defects in graphene created by irradiation or sample 
treatment [Tl fTTlfTQ] . 

In experiments, both symmetric and asymmetric vacancies are observed in unstressed graphene [7 , 
whereas in single layers of hexagonal Boron Nitride (hBN) only symmetric vacancies are observed [11 . 
Stable 5-8-5 divacancies are also observed |9]. The annihilation of the 5-7-7-5 SW defect in[5fb) (the 
heptagons share one side) 4 s after its creation is seen in Figures 3(c) and (d) of [7 . Our model predicts 
that SW under sufficient strain split in their two 5-7 pairs that move apart (cf Fig 6 of [13^), which 
has been observed very recently; cf Fig. 4(a) and (b) in [10]. 



4 Dislocation dipole pairs 

We study the evolution of configurations comprising two dislocation dipoles each in order to explain 
experimental observations by Meyer et al [7 . [6l^a) depicts an initial condition consisting of a 5-7-7-5 
SW defect adjacent to a rotated 7-5-5-7 defect (a, / = a/^/3^ h = 1/2 and H = 31/2 are the distances 
indicated in Figure [2]): 

U(x -x^,y- y^) - U(x - , ^ - % ) + U(x - x+, ^ - ^+) - U(x -x^.y- % ), (18) 
x^ = -0.3a, y^ = -0.7h + 2H, Xq = -0.3a - a, y^ = -0.7h + 2H, 

Xq = 0.3a + a, ^q" = 0.3h — ~ 0.15a, y^ = 0.5h. (19) 

Here U(x, y) is the edge dislocation ([16]) with origin of coordinates at a central atom of type A in Figure 
[Hand Burgers vector a (in units of the lattice constant a). \J{x^y) is an edge dislocation with Burgers 
vector b. To obtain U(x,^), we first consider the axes (x^y) rotated a 7r/3 angle from the axes ix^y). 
Next we form a 7-5-5-7 defect by combining a positive dislocation with Burgers vector (a, 0) centered at 
(xq",^q") and a negative dislocation with Burgers vector (— a,0) centered at (xq ,^^). Then the result 
is rewritten in the original coordinates (x^y). Now we add the 5-7-7-5 SW defect given by the first 
two terms in ([TS]) and complete that equation. The initial and boundary condition ([TS]) of the defect 
correspond to two dislocation dipoles having Burgers vectors along two different primitive directions 
and it is the same one as reported in Figures 3(h) and (i) of Meyer et al's experiments [7J. Under 
overdamped dynamics, this defect remains stable. As predicted in p^[T3]. the 7-5-5-7 defect is stable 
and this apparently stabilizes our pair of dislocation dipoles for the selected initial configuration. 
Other nearby configurations evolve to two octagons corresponding to a dipole comprising two edge 
dislocations with opposite Burgers vectors. As explained before, adsorbed atoms may be attached to 
the dangling bonds thereby eliminating these configurations and restoring the undisturbed hexagonal 
lattice. 

^ For an 80 keV incident electron, the maximum energy that can be transferred to a carbon atom is 15.8 
eV. This is below the threshold for knock-on damage (17 eV, corresponding to a beam energy of 86 keV) 
but sufficient to form multiple SW defects [7,22 . Migration of carbon atoms to empty neighboring sites has 
a significantly lower cost and therefore dislocation motion is much easier. The actual temperature rise in the 
suspended graphene specimen due to exposure by the electron beam in the microscope is minimal because 
the beam current density is small (^45 A/cm^) and the thermal conductivity of graphene is extremely high 
(>1000 W/mK). Thus the sample is not far from room temperature during the experiments [22j . 
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a) 



b) 



Fig. 6 (Color online) Defects corresponding to two dislocation dipoles: (a) Initial configuration of a pair of 
5-7-7-5 and 7-5-5-7 defects, (b) Final configuration. 




Fig. 7 (Color online) Metastable defect comprising three heptagons, three pentagons and one hexagon. It 
evolves to the defect-free lattice after 4 s. 
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The other hitherto unexplained defect configuration in Figures 3(j) and (k) of Ref. [7 is a metastable 
defect consisting of three pentagons and three heptagons. This defect appears with the fohowing initial 
and boundary condition depicted in Figure [71 

U(x -x-^.y- y-^) - \]{x - XQ,y -y^) ^ V{x -x^,y- y^) - V{x - Xq ,y - (20) 
Xq — — 0.3a + a, y^ = — 0.7/i, Xq = —0.3a, y^ = — 0.7/i, 

x+ = -0.3a - a, y^ = -0.7h - H, x^ = -0.3a, y^ = -OJh - H, (21) 

corresponding again to two dislocation dipoles all whose component dislocations have Burgers vectors 
directed along the x axis. Starting from a negative dislocation centered at {x^ ^y^) = (—0.3a, — 0.7/i), 
the first dipole adds a positive dislocation shifted one lattice constant to the left. The second dipole 
consists of a negative dislocation shifted vertically upwards a distance H = 31/2 = \^a/2 (1.5 times the 
hexagon side, or a/3/2 times the lattice constant) from {xq ^y^) and a positive dislocation which shifts 
horizontally to the right the previous one a distance equal to one lattice constant. Under overdamped 
dynamics, this defect disappears as the positive and negative dislocations comprising each dipole glide 
towards each other. See the movie in the Supplementary material. This agrees with Meyer et al's 
experimental observation 

5 Conclusions 

In summary, the proposed theory of defect dynamics in planar graphene regularizes continuum linear 
elasticity on a hexagonal lattice by replacing linear combinations of four appropriate difference op- 
erators acting on the displacement vector instead of partial derivatives thereof, and periodizes these 
operators along primitive directions. Far from defect cores where differences of the displacement vec- 
tor are sufficiently small, the resulting discrete equations seamlessly reduce to those of continuum 
linear elasticity. Adding large damping terms, these equations are solved with appropriate initial and 
boundary conditions consistent with the known solutions corresponding to edge dislocations and edge 
dislocation dipoles in linear elasticity. The numerical solutions of these equations explain the stability 
and evolution of several experimentally observed defects in suspended graphene sheets. Observed iso- 
lated defects are the cores of edge dislocations, dislocation dipoles or pairs of dipoles. Among them, 
isolated dislocations (pentagon- heptagon pairs), dislocation dipoles (symmetric vacancies, nonagon- 
pentagon pairs which are asymmetric vacancies, 5-8-5 divacancies, 5-7-7-5 Stone Wales defects and 
7-5-5-7 defects) and pairs of dislocation dipoles (a 5-7-7-5 SW defect adjacent to a 7-5-5-7 defect and a 
metastable defect comprising three pentagons, three heptagons and one hexagon). The theory correctly 
predicts stability or instability of these defects as observed in experiments. 
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